% =========================================================================
% bootstrap_tables.m
%  %Generates tables for Fulford and Schwartzman "The Benefits of
%  Commitment to a Currency Peg: The Gold Standard, National Banks, and the 1896 U.S. Presidential Election"
%
% Adapted from STICKY PRICES AND MONETARY POLICY: EVIDENCE FROM DISAGGREGATED U.S. DATA
% by Jean Boivin, Marc Giannoni, and Ilian Mihov 
% American Economic Review
% 
% 
% Inputs:
% data: Structure with data.y the aggregate indicator, data.x the
% information vector, data.year and data.month calendar variables.
% ts: Vector with identifying dates
% nf: number of factors
% norm: =1 if x is normalized so that all variances are equal to 1
% xlags: number of lags in factor extraction stage
% lags: number of lags in VAR
% Flags: number of lags in DFM
% bootstrap: = 1 if bootstrap is applied
% nrep1: number of repetitions in bootstrap for factor
% nrep2: number of repetitions in bootstrap for VAR for given factor draw.
% 
% =========================================================================


function out = bootstrap_tables(data,ts,nf,norm,xlags,lags,Flags,order,bootstrap,nrep1,nrep2)
out.tab = zeros(1,6);
out.Fst0 = [];
out.F1s = [];

if length(ts)<=nf %if not, error
%Housekeeping
 
x = data.x;
y = data.y(:,1);

year = data.year;
month = data.month;

y(isnan(y) & year>1900)=0;

iyy = find(1-isnan(y));
nfv = length(ts);
[T,~] = size(x);



if norm==1 %if norm=1, normalize x's by their std
    x = x./repmat(std(x),size(x,1),1);
end






%%%%%%%%%%%%%%%%%%% Estimate the observation equation %%%%%%%%%%%%%%%%%%%%%
%PC factors with all variables

[ Fst0, Fr0, lamr0, ebr0 ] = factor_st( x, nf, ts, xlags, Flags );


[ Fst0_bench, ~, ~, ~ ] = factor_st( x, 11, 80, 0 , 0);
    
if order ==1 
    fy0 = [y(iyy,1),Fst0(iyy,1)];
    py = 1;
    pf1 = 2;
    
elseif order ==2 
    
    fy0 = [Fst0(iyy,1),y(iyy,1)];
    py = 2;
    pf1 = 1;
 
end



%%%%%%%%%%%%%%%%%%%%%%%% Estimate VAR equation %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute lagged fy

lfy0 = ones(size(y(iyy,:),1),1);	% need 1's for irfbootstrap.m
for i = 1:lags
    lfy0 = [lfy0 lagn(fy0,i)];
end
fy0v  = fy0(lags+1:end,:);
lfy0 = lfy0(lags+1:end,:);

% VAR
b0   = olss(fy0v,lfy0);   
e0   = fy0v - lfy0*b0;

ty_gs = max(find(year==1900,2))-length(year)+length(fy0v);
b0_gs   = olss(fy0v(ty_gs:end,:),lfy0(ty_gs:end,:));   
Tg = length(fy0v(ty_gs:end,:));


k = size(fy0,2);
    
p = (size(lfy0,2)-1)/k;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%  Variance Decomposition and IRFs %%%%%%%%%%%%%%%%%%%

% Impulse response functions 

br  = zeros(size(b0,1),size(b0,2),nrep1*nrep2);
br_g  = zeros(size(b0,1),size(b0,2),nrep1*nrep2);


% Bootstrap for confidence intervals

if bootstrap ==1;
    
    
repetition = 0;
F1s = zeros(T-xlags,nrep1);
F1s_bench = zeros(T-xlags,nrep1);


for frep = 1:nrep1
    
    T = size(Fr0,1);
  	i = ones(T+1,1)+(T-1)*rand(T+1,1);
	i = round(i);
    
        
    xd = Fr0*lamr0' + ebr0(i(1:end-1),:);
    xd = [xd(1:xlags,:);xd];
    
    
[ Fst, ~, ~, ~ ] = factor_st( xd, nf, ts, xlags , Flags);

    F1s(:,frep) = Fst(:,1);
    
    
    


[ Fst_bench, ~, ~, ~ ] = factor_st( xd, 11, 80, 0 , 0 );


    F1s_bench(:,frep) = Fst_bench(xlags+1:end,1);    
    
    
if order ==1 
    fy = [y(iyy),Fst(iyy,1)];
    
elseif order ==2 
    
    fy = [Fst(iyy,1),y(iyy)];
 
end
    
    % compute lagged fy
    lfy = ones(size(fy,1),1);	% need 1's for irfbootstrap.m
    for i = 1:lags
        lfy = [lfy lagn(fy,i)];
    end
    fy  = fy(lags+1:end,:);
    lfy = lfy(lags+1:end,:);

    b = olss(fy,lfy);
    e = fy - lfy*b;
    
    T = size(fy,1);
    k = size(fy,2);
    
    
    b_g = olss(fy(T-Tg+1:end,:),lfy(T-Tg+1:end,:));
    e_g = fy(T-Tg+1:end,:)-lfy(T-Tg+1:end,:)*b_g;
    
    
    
    p = (size(lfy,2)-1)/k;

    
%    bxr = zeros(size(bx));

    btilda  = b;
    btilda_g = b_g;
%    bxtilda = bx;
           
    % Step 2a: IRF bootstrap using the bias-adjusted coeffs.
    for rep = 1:nrep2
        repetition = repetition+1;
        
        
        
%Full Sample        
    	i    = ones(T+1,1)+(T-1)*rand(T+1,1);
    	i    = round(i);
        
        
        lfyi = lfy(round(1+(T-1)*rand),:);
        
        
        
        lfyr = zeros(T+1,size(lfyi,2));
        fyr	 = zeros(T+1,size(e,2));
   
        for ii = 1:T+1
            fyi = lfyi*btilda + e(i(ii),:);
            fyr(ii,:) = fyi;
            lfyi = [1,fyi,lfyi(:,2:end-k)];
            lfyr(ii,:) = lfyi;
        end
        
        
        lfyr  = lfyr(1:T,:);
        fyr   = fyr(2:T+1,:);
        
        
%      	xr    = fyr*bxtilda + ex(i(1:end-1),:);
  %      bxr   = olss(xr,fyr);
    	br(:,:,repetition)    = olss(fyr,lfyr);
    	
      	
        
        %Only after gold standard

%Full Sample        
    	i    = ones(T+1,1)+(T-1)*rand(T+1,1);
    	i    = round(i);
        
        
        lfyi = lfy(round(1+(T-1)*rand),:);
        
        
        
        lfyr = zeros(T+1,size(lfyi,2));
        fyr	 = zeros(T+1,size(e,2));
   
        for ii = 1:T+1
            fyi = lfyi*btilda + e(i(ii),:);
            fyr(ii,:) = fyi;
            lfyi = [1,fyi,lfyi(:,2:end-k)];
            lfyr(ii,:) = lfyi;
        end
        
        
        lfyr  = lfyr(1:T,:);
        fyr   = fyr(2:T+1,:);
        
            	br(:,:,repetition)    = olss(fyr,lfyr);
    	
      	
%After gold standard        
        



        i    = ones(Tg+1,1)+(Tg-1)*rand(Tg+1,1);
    	i    = round(i);
        
        
        lfy_g = lfy(T-Tg+1:end,:);        
        lfyi = lfy_g(round(1+(Tg-1)*rand),:);
        
        
        
        lfyr = zeros(Tg+1,size(lfyi,2));
        fyr	 = zeros(Tg+1,size(e_g,2));
   
        for ii = 1:Tg+1
            fyi = lfyi*btilda_g + e_g(i(ii),:);
            fyr(ii,:) = fyi;
            lfyi = [1,fyi,lfyi(:,2:end-k)];
            lfyr(ii,:) = lfyi;
        end
        
        
        lfyr  = lfyr(1:Tg,:);
        fyr   = fyr(2:Tg+1,:);
        
            	br_g(:,:,repetition)    = olss(fyr,lfyr);
    	
      	
        
    end
    
    
corr = corrcoef(F1s(:,frep),F1s_bench(:,frep));
corrs(frep) = corr(1,2);
    
    frep
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Structural shocks
Sige = e0'*e0/(size(e0,1)-p*k-1);
smat = chol(Sige)';


%Calculate ``expectations shocks'' necessary to keep expectations constant:

fy1m = zeros(size(fy0));

fy1m(1:lags,:) = fy0(1:lags,:);
fy1m(1:lags,pf1) = 0;


fy1m_gs = zeros(size(fy0));

fy1m_gs(1:lags,:) = fy0(1:lags,:);
fy1m_gs(1:lags,pf1) = 0;

u0 = (inv(smat)*e0')';
u0_gs = u0;


for i = lags+1:length(fy0)
    
    u0(i-lags,pf1) = 1/smat(pf1,pf1)*(-[1 vec(flipud(fy1m(i-lags:i-1,:))')']*b0(:,pf1)-...
                                        (smat(pf1,:)*u0(i-lags,:)')'+smat(pf1,pf1)*u0(i-lags,pf1));
    
    fy1m(i,:) = [1 vec(flipud(fy1m(i-lags:i-1,:))')']*b0+(smat*u0(i-lags,:)')';
    
    
    u0_gs(i-lags,pf1) = 1/smat(pf1,pf1)*(-[1 vec(flipud(fy1m_gs(i-lags:i-1,:))')']*b0_gs(:,pf1)-...
                                        (smat(pf1,:)*u0_gs(i-lags,:)')'+smat(pf1,pf1)*u0_gs(i-lags,pf1));
    
    fy1m_gs(i,:) = [1 vec(flipud(fy1m_gs(i-lags:i-1,:))')']*b0_gs+(smat*u0_gs(i-lags,:)')';
    
end


if bootstrap ==1

    
fy1(1:lags,:) = fy0(1:lags,:);
fy1(1:lags,pf1) = 0;

fy1_gs(1:lags,:) = fy0(1:lags,:);
fy1_gs(1:lags,pf1) = 0;
    
    

    
for r=1:nrep1*nrep2
e0r = fy0v - lfy0*br(:,:,r);

e0r = e0;

Sige = e0r'*e0r/(size(e0r,1)-p*k-1);
smat = chol(Sige)';
u0 = (inv(smat)*e0r')';
u0_gs = u0;





for i = lags+1:length(fy0)
    
    
    u0(i-lags,pf1) = 1/smat(pf1,pf1)*(-[1 vec(flipud(fy1(i-lags:i-1,:))')']*br(:,pf1,r)-...
                                        (smat(pf1,:)*u0(i-lags,:)')'+smat(pf1,pf1)*u0(i-lags,pf1));
    
    fy1(i,:) = [1 vec(flipud(fy1(i-lags:i-1,:))')']*br(:,:,r)+(smat*u0(i-lags,:)')';
  
    
    u0_gs(i-lags,pf1) = 1/smat(pf1,pf1)*(-[1 vec(flipud(fy1_gs(i-lags:i-1,:))')']*br_g(:,pf1,r)-...
                                        (smat(pf1,:)*u0_gs(i-lags,:)')'+smat(pf1,pf1)*u0_gs(i-lags,pf1));
    
    fy1_gs(i,:) = [1 vec(flipud(fy1_gs(i-lags:i-1,:))')']*br_g(:,:,r)+(smat*u0_gs(i-lags,:)')';

    
    
    
    
    
end

fy1s(:,:,r) = fy1;
fy1s_gs(:,:,r) = fy1_gs;



    
    r
    end 
    

end



t_gs = max(find(year==1900,2))-xlags;

if bootstrap == 1;
out.tab = zeros(3,12);
ind = 2;

else
out.tab = zeros(1,12);
ind = 1;
end







%Relative std
out.tab(ind,1) = std(Fst0(t_gs:end,1))/std(Fst0(1:t_gs-1,1));




Fst0w = Fst0(:,1).*data.y(:,end)/data.y(80,end);

%Change around Sherman Purchase Act
out.tab(ind,2) = sum(Fst0w(find(year==1890 & month==7):find(year==1891 & month==7),1));

%Change after election
out.tab(ind,3) = sum(Fst0w(find(year==1896 & month==12):find(year==1897 & month==12),1));

%Change factor around 1893 panic
out.tab(ind,4) = Fst0w((year==1893 & month==7),1);

%Max abs change in fatctor after 1900
out.tab(ind,5) = max(abs(Fst0w(t_gs:end,1)));

yw = y(iyy,:).*data.y(iyy,end);
fy1mw = fy1m.*repmat(data.y(iyy,end),1,size(fy1m,2));
fy1m_gsw = fy1m_gs.*repmat(data.y(iyy,end),1,size(fy1m_gs,2));


out.y = fy1mw(:,py);

%Change around 1893 panic output
out.tab(ind,6) = exp(sum(yw(find(year==1893 & month==5)-length(year)+length(yw):find(year==1893 & month==10)-length(year)+length(yw))))-1;

%Change around 1893 panic, counterfactual
out.tab(ind,7) = exp(sum(fy1mw(find(year==1893 & month==5)-length(year)+length(fy1mw):find(year==1893 & month==10)-length(year)+length(fy1mw),py)))-1;

%Change around 1893 panic, Lucas counterfactual
out.tab(ind,8) = exp(sum(fy1m_gsw(find(year==1893 & month==5)-length(year)+length(fy1m_gsw):find(year==1893 & month==10)-length(year)+length(fy1m_gsw),py)))-1;


yy = y(iyy);

%Relative standard deviation before 1900
out.tab(ind,9) = std(fy1m(1:t_gs-1-length(year)+length(fy1m),py))/std(yy(1:t_gs-1-length(year)+length(yy)));

%Relative standard deviation after 1900 
out.tab(ind,10) = std(fy1m(t_gs-length(year)+length(fy1m):end,py))/std(yy(t_gs-length(year)+length(yy):end));

%Relative standard deviation before 1900, Lucas counterfactual
out.tab(ind,11) = std(fy1m_gs(1:t_gs-1-length(year)+length(fy1m_gs),py))/std(yy(1:t_gs-1-length(year)+length(yy)));

%Relateive standard deviation after 1900, Lucas counterfactual
out.tab(ind,12) = std(fy1m_gs(t_gs-length(year)+length(fy1m_gs):end,py))/std(yy(t_gs-length(year)+length(yy):end));

%Correlation with baseline
corr = corrcoef(Fst0(:,1),Fst0_bench(:,1));
out.tab(ind,13) = corr(1,2);


if bootstrap ==1
    
    l_ci_Fs = ceil(0.05*size(F1s,2));
    u_ci_Fs = floor(0.95*size(F1s,2));
    
    
    l_ci_y = ceil(.05*nrep1*nrep2);
    u_ci_y = floor(.95*nrep1*nrep2);
   
   sdF1s_rat = std(F1s(t_gs:end,:))./std(F1s(1:t_gs-1,:));
   sdF1s_rat = sort(sdF1s_rat);
   out.tab(1,1) = sdF1s_rat(l_ci_Fs);
   out.tab(3,1) = sdF1s_rat(u_ci_Fs);
   
   F1sw = F1s.*repmat(data.y(:,end),1,size(F1s,2))/data.y(80,end);

   sums = sum(F1sw(find(year==1890 & month==7):find(year==1891 & month==7),:),1);  
   sums_sort = sort(sums,2);
   out.tab(1,2) = sums_sort(l_ci_Fs);
   out.tab(3,2) = sums_sort(u_ci_Fs);
   
   
   
   sums = sum(F1sw(find(year==1896 & month==12):find(year==1897 & month==12),:),1);  
   sums_sort = sort(sums,2);
   out.tab(1,3) = sums_sort(l_ci_Fs);
   out.tab(3,3) = sums_sort(u_ci_Fs);
   
      
   
   F1s_sort = sort(F1sw,2);   
   out.tab(1,4) = F1s_sort((year==1893 & month==7),l_ci_Fs);
   out.tab(3,4) = F1s_sort((year==1893 & month==7),u_ci_Fs);
   
   max_sort = sort(max(abs(F1sw(t_gs:end,:)),[],1));
   out.tab(1,5) = max_sort(l_ci_Fs);
   out.tab(3,5) = max_sort(u_ci_Fs);      
   
   out.tab([1,3],6) = out.tab(2,6);
   
   
   fy1sw = fy1s.*repmat(data.y(iyy,end),[1,size(fy1s,2),size(fy1s,3)]);
   
   
   
   sums =  sum(fy1sw(find(year==1893 & month==5)-length(year)+length(fy1m):find(year==1893 & month==10)-length(year)+length(fy1m),py,:),1);
   sums = sort(squeeze(sums));
   out.tab(1,7) = exp(sums(l_ci_y))-1;
   out.tab(3,7) = exp(sums(u_ci_y))-1;
   
   
   fy1s_gsw = fy1s_gs.*repmat(data.y(iyy,end),[1,size(fy1s,2),size(fy1s,3)]);
   
   sums_g =  sum(fy1s_gsw(find(year==1893 & month==5)-length(year)+length(fy1m):find(year==1893 & month==10)-length(year)+length(fy1m),py,:),1);
   sums_g = sort(squeeze(sums_g));
   out.tab(1,8) = exp(sums_g(l_ci_y))-1;
   out.tab(3,8) = exp(sums_g(u_ci_y))-1;

    
   sdy1s_1 = std(fy1s(1:t_gs-1-length(year)+size(fy1s,1),py,:),[],1);
   sdy1s_1 = sort(squeeze(sdy1s_1));   
   out.tab(1,9) = sdy1s_1(l_ci_y)/std(yy(1:t_gs-1-length(year)+length(yy)));
   out.tab(3,9) = sdy1s_1(u_ci_y)/std(yy(1:t_gs-1-length(year)+length(yy)));
   
   
   sdy1s_2 = std(fy1s(t_gs-length(year)+size(fy1s,1):end,py,:),[],1);
   sdy1s_2 = sort(squeeze(sdy1s_2));   
   out.tab(1,10) = sdy1s_2(l_ci_y)/std(yy(t_gs-length(year)+length(yy):end));
   out.tab(3,10) = sdy1s_2(u_ci_y)/std(yy(t_gs-length(year)+length(yy):end));         
   
   
   
   sdy1s_1 = std(fy1s_gs(1:t_gs-1-length(year)+size(fy1s,1),py,:),[],1);
   sdy1s_1 = sort(squeeze(sdy1s_1));   
   out.tab(1,11) = sdy1s_1(l_ci_y)/std(yy(1:t_gs-1-length(year)+length(yy)));
   out.tab(3,11) = sdy1s_1(u_ci_y)/std(yy(1:t_gs-1-length(year)+length(yy)));
   
   
   sdy1s_2 = std(fy1s_gs(t_gs-length(year)+size(fy1s,1):end,py,:),[],1);
   sdy1s_2 = sort(squeeze(sdy1s_2));   
   out.tab(1,12) = sdy1s_2(l_ci_y)/std(yy(t_gs-length(year)+length(yy):end));
   out.tab(3,12) = sdy1s_2(u_ci_y)/std(yy(t_gs-length(year)+length(yy):end));         
   
   corrs = sort(corrs);
   
   out.tab(1,13) = corrs(l_ci_Fs);
   out.tab(3,13) = corrs(u_ci_Fs);
   
   
   
   
   out.F1s = F1s;
    out.fy1s = fy1s;
    out.fy1sw = fy1sw;
   
end


   out.Fst0 = [zeros(size(x,1)-size(Fst0,1),size(Fst0,2));Fst0];


end


